ACID RAIN: MICROPHYSICAL MODEL 


By A. N. Dingle* 


BASIC EQUATIONS 


When HCl and H 2 O vapors are present together, they jointly determine the 
equilibrium vapor pressure of each component over an HCl^- droplet. It is 
assumed that the two vapors may be treated as if they diffuse independently. 
Thus, the change of mass of a solution droplet may be expressed as 


dm = dm- + dm^ 
r 1 2 


( 1 ) 


where the subscripts indicate the respective components; i.e., H 2 O and HCl. 
Then, analogously with equation (1), 
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The energy balance of the droplet may be expressed by 
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where Q^j;, the internal energy increase, is defined as 


” ®r^r dt 


Ql, the latent energy release, is defined as 






Ql^ = L2(diU2/dt 


Qk, the conductive heat transfer from drop, is defined as 
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Q^, the radiative transfer from drop, is defined as 
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Q^, the surface increase energy, is defined as 
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Qj 4, mixing, is defined as 
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Qp, frictional conversion, is defined as 
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and Qj), the heat of dilution, is defined as 
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A detailed discussion of the energy terms is given in reference 1 and 
need not be repeated here. A study of the relative magnitudes of these terms 
for different sizes of droplets of 4.0 molal HCl^q in a 1-m/sec updraft (table 
I) shows that only Qp is definitely negligible in the cloud droplet size 
range. This allows simplification of equation (4) to 

'’t ■ \ + '’k + 'Jr - ‘Ir + 


By means of the definitive expressions, the droplet temperature elevation may 
be written as 
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In addition, for each time step, the mass conservation accounts are 
kept. Assuming that there is no dilution of the air parcel by turbulent 
mixing with environmental air. 
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The environmental temperature and pressure changes from one time step to an- 
other are given by 
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and 
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where the summation in the denominator represents the heat capacity of the 
liquid carried with the rising air parcel. Equations (2), (3), and (6) for 
each droplet size category and equations (7), (8), (9), and (10) constitute 
the explicit model for two volatile components. 
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NUMERICAL INTEGRATION TECHNIQUES 


The numerical technique used for the integration of the microphysical 
equations is Hamming's modified predictor-corrector method started by a 
Runge-Kutta procedure. This method has been described in reference 1 and is 
given in reference 2. At present, there are 2n + 2 first-order ordinary 
differential equations to be solved simultaneously, where n is the number 
of droplet size classes. In terms of the variables used in reference 2, 
yj through y^ are the derivatives with respect to time of the mass of 
water in a droplet in a specific size class: 
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y^+2 through y^^ are the derivatives with respect to time of the 
HCl mass of a droplet in a specific size class: 


2n 

n+1 


diUo 


dt 


47rrD^F V /e 

^ ''2 / *2 


69 



is the derivative with respect to time of cloud air temperature: 




dT 

^ _ a 

^2n+l ” dt 



( 11 ) 


where n£ in the denominator on the right-hand side of equation (11) refers 

not to the number of size classes but to the number of droplets in size 

class i and where y' ^ is the derivative with respect to time of the 
, , . -^2n+2 

cloud air pressure; 
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The droplet distribution is initialized by assuming a starting pressure 
or height (y2n+2^®^^5 temperature (y 2 n+ 1 ^ 0 )), relative humidity, and vapor 
pressure of HCl of the center of the cloud. The small droplets are assumed 
to be in equilibrium with the environment and the large droplets near equilib- 
rium. Under these assumptions, the mass of H 2 O (yi(0) and mass of 

HCl (yxi+l(0) Y2n^^^^ each droplet are determined. Integration of the 

microphysical equations now begins. 

One of the benefits of the Hamming method is the fact that the local trun- 
cation error (LTE) is calculated at each increment of integration or time step 
and can be used to modify this time step. If the sum of the LTE*s for each 
equation exceeds a given tolerance limit specified by the user, a new time 
step that is one-half the length of the previous time step is selected and the 
integration is restarted at the last accepted integrated values. This halving 
of the time step continues until the new values are within the error limit or 
until the number of halvings exceeds a specified limit. If the sura of LTE*s 
is between the tolerance limit and one-fiftieth of the given tolerance limit, 
the calculated values are accepted and integration continues. Finally, if the 
Stan of the LTE’s is less than one-fiftieth of the given tolerance limit, the 
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calculated values are accepted, the time step is doubled, and the integration 
continues. With this time step adjustment, the integration is more accurate 
and efficient than if a fixed time step were used. 

A computational instability was encountered during the integration of the 
microphysical equations. Some droplets were found to fluctuate wildly about 
their equilibrium values. The vapor pressures of HCl and H2O in the environ- 
ment may not equal the equilibrium vapor pressures at the droplet surface 
which are a function of temperature, molality, and radius. The existence of 
a vapor gradient allows the molality of HCl or H2O (or both) independently of 
each other. The consequent condensation/ evaporation of vapors may overcompen- 
sate the change in molality necessary for the droplet to reach equilibrium with 
the environment. The overcompensation of molality causes the vapor gradient 
to change sign, and during the next step, growth is in the opposite direction. 
To correct for the resultant instability or oscillation of the droplets about 
their equilibrium, the following scheme was developed. 

Consider, for the moment, one diffusing vapor. The vapor difference at 
the droplet surface is some value B: 
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If the droplet goes to equilibrium, B goes to zero. The change of vapor 
difference for a droplet that goes to equilibrium during the time step is 
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Grouping terms, one obtains 
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( 12 ) 


The vapor gradient is not allowed to collapse far beyond zero in a time step. 
Therefore, the change in vapor pressure Ae defines the maximum allowable 
change in vapor pressure. The maximum change in vapor pressure can then be 
used to find the maximum allowable change in droplet molality during a time 
step. The derivation of the maximum change in molality follows. 
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The vapor pressure to a droplet can be expressed as a function of drop- 
let molality, temperature, and radius. 


e = f(y,T,r) 


Differentiating, one finds that 
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Having the equations for the vapor pressure of each component over a flat 
solution as a function of T and y (ref. 1), it is an easy matter to 
differentiate and find that 




and 


9T 


Now, for a droplet of radius r, the vapor pressure at the surface is 
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Then, 
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Transforming dr dm, one obtains 
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Substituting equations (14a), (14b), and (14c) into equation (13) yields 


de = A dy + F dT + G dm 
and in finite difference form, 

Ae = A Ay + F AT -h G Am 
de 

= B + ^ AT (15) 


From the definition of molality, one finds that 


Ay 




where m^ is the solvent and m 2 is the solute. 
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The relative growth rate (dm 2 /dm 2 ^) can be determined from the growth 
equations (2) and (3): 
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Then, 


Ay = y 




(16) 


Substituting equation (16) into equation (15) and solving for Ay> 
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where d(e 3 )/dt and dT^/dt are estimated by backward differencing. Then, 


F AT 



The maximum allowable change of molality determined from equation (17), 
Ayi, applies to the solvent. Similar calculations for the solute result in 
another value, Ay 2* With the smaller value of the two, the maximum allow- 
able values of A^l A^2 can be found in equation (16). This method is 

approximate. 
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However, one would not want to prevent the droplet from crossing the 
equilibrium line when the natural process is to do so. If the growth of the 
droplets becomes unstable, it will fluctuate wildly across the equilibrium 
line. Therefore, a domain about the equilibrium line is defined. If the 
droplet tends to grow and/or evaporate quickly outside that domain, it is 
assumed that there is a numerical instability. For the current computation, 
if the change of mass of either component computed by the numerical integra- 
tion is less than three times the maximum allowable change estimated above, 
the integrated value is assumed correct. In practice, the factor of 3 is 
quite arbitrary; however, it is found that the ratio of calculated Am to 
the estimated maximum allowable A^i is either «1 or »10. If the droplet 
growth is constrained in this manner, after the first few seconds on integra- 
tion, only the two smallest size classes need adjustment and the model is 
stable. 


RESULTS 


To assess the situation in which droplet growth should be most favored, 
the microphysical model has been used to simulate the case of a ground cloud 
(GC) without dilution by entrainment and without precipitation. The updraft 
speeds used (table II) were chosen to approximate the observed behavior of 
several Titan ground clouds (ref. 3) as shown in figure 1. Thus, in this 
case, the GC reaches a height of 1000 meters at t - 250 seconds and stops 
rising after t = 500 seconds at a level of 1500 meters. 

The initial GC temperature (at t = 90 seconds) is set at 25*^ C, and the 
initial HCl vapor pressure is 3.8 dyn/cm^, corresponding to a concentration 
of 4 p/m at the 950-millibar level. The dry particle distribution was taken 
from the suggestion of Lala (ref. 4). It is a bimodal distribution of the 
following form. 
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The various parameters suggested by Lala (ref. 4) are as follows. 


^2 " N ^/8000 

Y] N = 1.5 X 10^ particles/om'^ 
r=0.01 

a T = a ^ 
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= 2,0 micrometers 

r - - 0.05 micrometer 
gl 


r ^ ~ 1.0 micrometer 

g2 

The distribution is plotted in figure 2. 

The course of the relative humidity H for the period 90 seconds < t 
< 720 seconds is shown in figure 3 for initial humidity values of 80, 90, and 
100 percent. The effect of rapid initial droplet growth in the 100-percent 
case is to reduce the initial humidity quickly to a minimum of 98 percent 
at t = 99 seconds and then, with continued rising and cooling of the cloud, 
to approach the nearly constant value of 99.7 percent after 200 seconds. 

The as 3 rmptotic humidity value decreases by about 0.1-percent steps to the 
90- and 80-percent cases. 

Figure 4 shows the course of the HCl vapor pressure. It is reduced to 
lO""^ dyn/cm^ at about t = 300 seconds for the 100-percent case, about 100 
seconds later for the 90-percent case, and not until t 500 seconds for the 
80-percent case. The droplet molalities decrease to 0.1 at t = 225, 300, 
and 500 seconds for the 100- , 90- , and 80-percent cases, respectively 
(fig. 5). 


The droplet-size spectra after 720 seconds are shown in figure 6. The 
outstanding feature of this set of results is the enhancement of the growth 
of the largest droplets as the initial humidity is lowered. The curves are 
truncated at 7 micrometers, although the largest drops formed exceed 20 mi- 
crometers in radius. These are, however, very few in number. The bifurca- 
tion in the 100-percent relative humidity size distribution curve is an 
artifact produced in this case by the mechanics of calculating dN/dr. The 
equivalent cloud liquid water content figures at 720 seconds are, respective- 
ly, 2.8, 2.1, and 1.6 g/m^ with mean droplet radii of 3.8, 3.75, and 3.4 
micrometers. 
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An additional result of this study is that the droplet molalities are 
nearly constant across the size distribution at each value of t, a fact 
which will be helpful in parameterizing the microphysical processes for the 
submesoscale dynamic model. A summary of results at t = 720 seconds is 
given in table III. 


REFERENCES 


1. Dingle, A. N.: Rain Scavenging of Solid Rocket Exhaust Clouds. NASA 

CR-2928, 1978. 

2. System/360 Scientific Subroutine Package (360A-CM-03X) Version 3 Program- 
mer's Manual. International Business Machines Corp., Technical Publica- 
tions Dept. H20-0205-3 (White Plains, N.Y.), 1970, p. 545. 

3. Stephens, J. B.; and Stewart, R.: Rocket Exhaust Effluent Modeling for 

Tropospheric Air Quality and Environmental Assessments; Preliminary Re- 
port for NASA Space Shuttle Tropospheric Environmental Effects Meeting, 
Langley Research Center (Hampton, Va. 24026), Feb. 1975. 

4. Lala, G. G.; First draft for Position Paper on the Potential of Inadver- 
tent Weather Modification of the Florida Peninsula Resulting From Neu- 
tralization of Space Shuttle Solid Rocket Booster Exhaust Clouds, by 

E. Ballay, L. Bosart, E. Droessler, J. Jiusto, G. G. Lala, V. Mohnen, 

V. Schaefer, and P. Squires. Prepared for NASA Langley Research Center 
by Institute of Man and Science (Rensselaerville, N.Y.), Feb. 1978. 


77 



TABLE I.- ESTIMATES OF THE ENERGY TERMS 


[in relation to Q-j for the HCI/H 2 O system at p “ 


Energy term 

r = 0.1 ym 

r = 1. 

0 ym 

r == 10 

Ql/Qt ^*^Qr/Qt^ 

1.46 X 10^ 

4.28 

X 10^ 

1.95 X lo3 

Qr/Qt (e = 1 ) 

.92 

9.7 


3.7 

Qr/Qi 

8.9 

2.6 


.012 

Qy/QT 

.006 

1.05 


.24 

Qp/Qx 

8 X 10~15 

1.2 X 

10-12 

1.2 X 10"10 

Qp/Qi 

94 

276 


13 


TABLE II.- UPDRAFT SPEEDS AND HEIGHTS FOR 
STABILIZED GROUND CLOUD 


Time, 

sec 

Updraft speed, 
m/sec 

Height, m 

90 to 

250 

4 

360 to 1000 

250 to 

500 

2 

1000 to 1500 

500 to 

720 

0 

1500 
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TABLE III.- VARIABLES AT TIME = 720 SECONDS 


Variable 

Initial 

80 

humidity, 

90 

percent 

100 

Temperature, K 

290.99 

292.45 

293.73 

Relative humidity, percent . . . 

99.59 

99.65 

99.70 

HCl vapor pressure, dyn/cm^ , . , 

0.001 

0.001 

0.001 

Droplet molality 

0.09 

0.07 

0.05 

Liquid water content, g/m^ • . . 

1.56 

2.14 

2.80 

Mean radius, ym 

3.40 

3.74 

3.79 

Standard deviation, ym 

0.45 

0.62 

0.31 

Largest droplet, ^m 

19.56 

20.91 

17.04 
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Figure 1.- Cloud rise and stabilization heights in eight Titan III launches. 



PARTICLE RADIUS, Mm 

Figure 2.- Particle size distribution of AI 2 O 3 in the ground cloud (after 

ref. 4). 
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I STARTING RH: 100 PERCENT 



Figure 3.- Course of the relative humidity from t = 90 seconds to t = 720 
seconds for initial humidity values of 80, 90, and 100 percent. 



Figure 4.- Course of the HCl vapor pressure. 
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DROPLET MOLALITY 



0 100 200 300 400 500 600 700 800 

TIME, SEC 


Figure 5*- Change of droplet molality with time. 



RADIUS, /urn 

Figure 6.- Droplet-size spectra after 720 seconds (truncated at r = 7 

micrometers) . 
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